package org.javanum.number;

import java.math.BigDecimal;
import java.math.MathContext;

/**
 *
 * <p>This class consists exclusively of static utility methods that
 * operate on BigDecimal and BigInteger classes.</p>
 *
 * <p>All methods in this class are allowed to throw a
 * NullPointerException if null is passed as an argument.</p>
 *
 * Author: Scott Fines
 * Date: Dec 11, 2009
 *
 */
public class BigMath {

    private BigMath(){}

    public static volatile int MAX_ITERATIONS = 50;
    public static final BigDecimal TWO = new BigDecimal("2");
    public static final BigDecimal oneHalf = BigDecimal.ONE.divide(TWO);


    /**
     * <p>Approximates the square root of the specified BigDecimal
     * to the accuracy defined by the specified MathContext. </p>
     *
     * <p>This method uses an <i>iterative algorithm</i>to approximate
     * the root. Iterative algorithms perform a series of calculations
     * repeatedly until the desired accuracy is reached.
     * In some situations, the desired accuracy may never be reached,
     * which can result in infinite loops. To avoid this situation,
     * this method uses the {@link #MAX_ITERATIONS} field to determine
     * a maximum number of steps to be taken before the algorithm is
     * forcibly terminated. In this situation, an ArithmeticException
     * will be thrown.</p>
     *
     * <p>If the user is confident that the approximation will
     * in fact terminate, then setting {@link #MAX_ITERATIONS}
     * to a higher value will allow greater accuracy. However, setting this
     * value too high may result in poor overall performance.</p>
     * 
     * @param value the BigDecimal whose square root should be evaluated
     * @param rootContext the MathContext specifying the desired accuracy
     * @return an approximation of the square root of value,
     *          to the accuracy specified by rootContext
     * @throws ArithmeticException if the number of steps required
     *          to converge is greater than {@link #MAX_ITERATIONS},
     *          or if value.signum()<0
     */
    public static BigDecimal bigRoot(BigDecimal value,
                                     MathContext rootContext){
        //TODO don't use 2*precision for every calculation!
        if(value.signum()<0)
            throw new ArithmeticException(
                    "Cannot compute the square root of a negative number");
        if(equalToWithinPrecision(value,BigDecimal.ZERO,
                                        rootContext.getPrecision()))
            return BigDecimal.ZERO;
        MathContext divideContext =
                new MathContext(rootContext.getPrecision()*2);

        int stepCount=0;
        BigDecimal tolerance =
                new BigDecimal("1").
                        movePointLeft(rootContext.getPrecision());
        BigDecimal p0 = new BigDecimal(Math.sqrt(value.floatValue()));
        if(BigDecimal.ZERO.equals(p0))
            p0 = value;
        BigDecimal p;

        do{
            BigDecimal t1 = p0.divide(TWO,divideContext);
            BigDecimal t2 = TWO.multiply(p0);
            BigDecimal v1 = value.divide(t2,divideContext);
            p = t1.add(v1);
            if(p.subtract(p0).abs().compareTo(tolerance)<0){
                return p; 
            }
            p0=p;
        }while(stepCount< MAX_ITERATIONS);

        throw new ArithmeticException(
                "The desired accuracy could not be attained within the " +
                "maximum number of allowed steps.");
    }

    /**
     * <p>Determines whether two BigDecimals are equal to within a
     * specified precision.</p>
     *
     * @param value1 a number to compare
     * @param value2 the other number to compare
     * @param precision the number of decimal places of precision
     *          to compare the two numbers with
     * @return true if two numbers are identical in the first
     *          {@code precision} number of decimal places
     */
    public static boolean equalToWithinPrecision
            (BigDecimal value1, BigDecimal value2, int precision){
        BigDecimal tol = BigDecimal.ONE.movePointLeft(precision);
        return value1.subtract(value2).abs().compareTo(tol)<0;
    }


    public static BigDecimal sin(BigDecimal input,MathContext precision){
        return null;
    }

    public static BigDecimal cos(BigDecimal input,MathContext precision){
        return null;
    }

    public static BigDecimal tan(BigDecimal input,MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal sec(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal csc(BigDecimal input,MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal cot(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal asin(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal acos(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal atan(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal asec(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal acsc(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal acot(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal sinh(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal cosh(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal tanh(BigDecimal input,MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal asinh(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal acosh(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal atanh(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal asech(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal acsch(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal acoth(BigDecimal input, MathContext precision){
        //TODO implement!
        return null;
    }

    public static BigDecimal exp(BigDecimal input, MathContext precision){
        BigDecimal exp = new BigDecimal(Math.exp(input.doubleValue()));
        BigDecimal oldExp;
        int interPrec = precision.getPrecision();
        BigDecimal tol = BigDecimal.ONE.movePointLeft(interPrec);
        interPrec+=(int)Math.log(interPrec);
        MathContext interContext = new MathContext(
                            interPrec,precision.getRoundingMode());
        BigDecimal newX = input.add(BigDecimal.ONE);
        do{
            oldExp = exp;
            exp = newX.subtract(ln(oldExp,interContext)).multiply(oldExp);
        }while(exp.subtract(oldExp).abs().compareTo(tol)>=0);

        return exp.round(incrementPrecision(precision));
    }

    private static MathContext incrementPrecision(MathContext precision){
        return new MathContext(
                precision.getPrecision()+1,precision.getRoundingMode());
    }

    private static MathContext incrementPrecision(MathContext precision,int increment){
        return new MathContext(
                precision.getPrecision()+increment,precision.getRoundingMode());
    }


    /**
     * <p>Returns the number which is PI, approximated to the
     * specified number of decimal places.</p>
     *
     * <p>Note that this method will always return pi approximated
     * to <i>at least</i> the precision requested, but may return
     * an approximation which is more precise.</p>
     *
     * <p>Note that this implementation caches the previous
     * results, and compares the precision of that previous
     * result with the requested precision. If the previous
     * result is sufficiently precise, then it is returned without
     * calculation. Otherwise, pi is calculated to the correct
     * precision and stored for later use. Therefore, it is
     * possible that the returned value is much more precise than
     * what is requested.</p>
     *
     * <p>This method is thread-safe.</p>
     *
     * <p>This implementation is thread-safe, but some of the
     * optimizations taken may not be maintained across threads.
     * That is, two separate threads may perform the same calculation,
     * rather than making use of the caching mechanism.</p>
     *
     * @param precision the desired precision
     * @return the value of pi which is approximated to {@code precision}
     *          decimal places.
     */
    public static BigDecimal pi(MathContext precision){
        return Pi.instance.pi(precision);

    }

    private enum Pi{
        instance;
        private volatile BigDecimal pi = BigDecimal.valueOf(Math.PI);
        private volatile MathContext calcPrecision = new MathContext(17);
        private volatile BigDecimal an= BigDecimal.ONE,
                                    bn,
                                    tn =BigDecimal.ONE.
                                        divide(new BigDecimal("4.0")),
                                    pn = BigDecimal.ONE;

        public BigDecimal pi(MathContext precision){
            if(calcPrecision.getPrecision()>=precision.getPrecision())
                return pi.round(precision);

            return calculateAndStorePi(precision.getPrecision()+1);
        }

        private synchronized BigDecimal calculateAndStorePi(int precision) {
            //TODO don't use 2*precision for every calculation!
            MathContext prec = new MathContext(2*precision);
            if(bn==null){
                BigDecimal sqrtTwo =
                        bigRoot(TWO,prec);
                bn = BigDecimal.ONE.divide(sqrtTwo,prec);
            }
            BigDecimal tol = BigDecimal.ONE.movePointLeft(precision);
            do{
                BigDecimal a = an;
                an = an.add(bn).multiply(oneHalf);
                bn = bigRoot(a.multiply(bn),prec);
                tn = tn.subtract(a.subtract(an).pow(2).multiply(pn));
                pn = TWO.multiply(pn);
            }while(an.subtract(bn).abs().compareTo(tol)>0);

            calcPrecision = new MathContext(precision);
            pi = an.add(bn).pow(2).divide(
                    new BigDecimal("4.0").multiply(tn),calcPrecision);

            return pi;
        }
    }

    /**
     * <p>Returns the number which is ln(2), approximated to the
     * specified number of decimal places.</p>
     *
     * <p>ln(2) is defined as the number such that {@code e^ln(2)=2},
     * where {@code e=2.71828...} is the base
     * of the natural logarithm</p>
     *
     * <p>Note that this method will always return ln(2) approximated
     * to <i>at least</i> the precision requested, but may return
     * an approximation of higher precision.</p>
     *
     * <p>Note that this implementation caches the previous
     * results, and compares the precision of that previous
     * result with the requested precision. If the previous
     * result is sufficiently precise, then it is returned without
     * calculation. Otherwise, ln(2) is calculated to the correct
     * precision and stored for later use. Therefore, it is
     * possible that the returned value is much more precise than
     * what is requested. To ensure that the returned number is
     * always exactly the precision asked for,
     * call {@link #ln2(java.math.MathContext, boolean)} instead.</p>
     *
     * <p>This method is thread-safe.</p>
     *
     * <p>This implementation is thread-safe, but some of the
     * optimizations taken may not be maintained across threads.
     * That is, two separate threads may perform the same calculation,
     * rather than making use of the caching mechanism. However,
     * the caching mechanism satisfies a <i>happens before</i>
     * relationship.</p>
     *
     * @param precision the desired precision
     * @return the value of ln(2) which is
     *          approximated to <i>at least</i>
     *          {@code precision} decimal places.
     */
    public static BigDecimal ln2(MathContext precision){
        return ln2(precision,false);
    }

    /**
     * <p>Returns the number which is ln(2), approximated to the
     * specified number of decimal places.</p>
     *
     * <p>ln(2) is defined as the number such that {@code e^ln(2)=2},
     * where {@code e=2.71828...} is the base
     * of the natural logarithm</p>
     *
     * <p>Note that this method will always return ln(2) approximated
     * to <i>at least</i> the precision requested, but may return
     * an approximation of higher precision.</p>
     *
     * <p>Note that this implementation caches the previous
     * results, and compares the precision of that previous
     * result with the requested precision. If the previous
     * result is sufficiently precise, then it is returned without
     * calculation, rounding as necessary to return a number which has
     * exactly the precision specified.
     * Otherwise, ln(2) is calculated to the correct
     * precision and stored for later use. To avoid the extra
     * step of rounding, call {@link #ln2(java.math.MathContext)}
     * instead.</p>
     *
     * <p>This method is thread-safe.</p>
     *
     * <p>This implementation is thread-safe, but some of the
     * optimizations taken may not be maintained across threads.
     * That is, two separate threads may perform the same calculation,
     * rather than making use of the caching mechanism. However,
     * the caching mechanism satisfies a <i>happens before</i>
     * relationship.</p>
     *
     * @param precision the desired precision
     * @return the value of ln(2) which is approximated to
     *          <i>exactly</i> {@code precision} decimal places.
     */
    public static BigDecimal ln2(
            MathContext precision, boolean shouldRound){
        if(shouldRound)
            return ln2.instance.compute(precision).round(precision);
        
        return ln2.instance.compute(precision);
    }

    private static enum ln2{
        instance;
        private volatile BigDecimal ln2 = new BigDecimal(Math.log(2));
        private volatile MathContext currentContext = new MathContext(17);

        public BigDecimal compute(MathContext precision){
            if(currentContext.getPrecision()>=precision.getPrecision())
                return ln2;
            else
                return calculate(precision);
        }

        private synchronized BigDecimal calculate(MathContext precision) {
            MathContext prec = new MathContext(precision.getPrecision()+1,
                                               precision.getRoundingMode());
            int m = getM(prec.getPrecision());
            BigDecimal twoUp = oneHalf.pow(m-2);
            BigDecimal agm = AGM.getInstance(BigDecimal.ONE,twoUp).compute(prec);
            BigDecimal bottom = new BigDecimal(m).multiply(TWO).multiply(agm);
            BigDecimal pi = pi(prec);
            ln2 = pi.divide(bottom,prec);

            currentContext = new MathContext(precision.getPrecision(),
                                            precision.getRoundingMode());
            return ln2;
        }

        private int getM(int precision) {
            return 2*precision+2;
        }

    }

    public static BigDecimal ln(BigDecimal x, MathContext precision){
        if(BigDecimal.ZERO.compareTo(x)>0)
            throw new ArithmeticException("Cannot compute the natural logarithm" +
                    "of a negative number");
        if(BigDecimal.ZERO.compareTo(x)==0)
            throw new ArithmeticException("ln(0)=-infinity cannot be computed");
        if(BigDecimal.ONE.equals(x))
            return BigDecimal.ZERO;

        int m = 2*precision.getPrecision()+2;
        BigDecimal twoUp = oneHalf.pow(m-2);
        int interimPrecision = precision.getPrecision();
        interimPrecision+=(int)Math.log(interimPrecision);
        MathContext prec = new MathContext(interimPrecision,
                                    precision.getRoundingMode());
        boolean lessThanOne=false;
        if(x.compareTo(BigDecimal.ONE)<0){
            lessThanOne=true;
            BigDecimal recip = BigDecimal.ONE.divide(x,prec);
            twoUp = twoUp.divide(recip,prec);
        }else{
            twoUp = twoUp.divide(x,prec);
        }

        BigDecimal agm = AGM.getInstance(BigDecimal.ONE,twoUp).compute(prec);
        BigDecimal bot = TWO.multiply(agm);
        BigDecimal main = pi(prec).divide(bot,prec).
                            subtract(new BigDecimal(m).
                            multiply(ln2(prec))).
                            round(incrementPrecision(precision));
        if(lessThanOne)
            return main.negate();
        else
            return main;
    }

    public static BigDecimal e(MathContext precision) {
        return E.instance.e(precision);
    }

    private static enum E {
        instance;

        private volatile BigDecimal e = new BigDecimal(Math.E);
        private volatile MathContext currentContext = new MathContext(17);
        private volatile BigDecimal p1 = e;
        public BigDecimal e(MathContext precision){
            if(currentContext.getPrecision()>=precision.getPrecision())
                return e;

            return calculateAndStoreE(precision);
        }

        private synchronized BigDecimal calculateAndStoreE(MathContext precision) {
            BigDecimal p0;
            BigDecimal tol =
                    BigDecimal.ONE.movePointLeft(precision.getPrecision());
            int interPrec = precision.getPrecision();
            interPrec+=(int)Math.log(interPrec);
            MathContext prec = new MathContext(
                    interPrec,precision.getRoundingMode());

            do{
                p0=p1;
                p1 = TWO.subtract(ln(p0,prec)).multiply(p0);
            }while(p1.subtract(p0).abs().compareTo(tol)>=0);

            currentContext = precision;
            e = p1;

            return e.round(incrementPrecision(precision));
        }

    }

    private static class AGM{
        private volatile BigDecimal a1;
        private volatile BigDecimal b1;
        private volatile MathContext context;

        private AGM(BigDecimal a0, BigDecimal b0){
            this.a1=a0;
            this.b1=b0;
            this.context = new MathContext(0);
        }
        
        public static AGM getInstance(BigDecimal a0, BigDecimal b0){
            return new AGM(a0,b0);
        }
        public BigDecimal compute(MathContext precision){
            if(precision.getPrecision()<=context.getPrecision())
                return a1;
            else
                return computeMean(precision);
        }

        private synchronized BigDecimal computeMean(MathContext precision) {
            BigDecimal a0;
            BigDecimal b0;
            BigDecimal TOL = BigDecimal.ONE.movePointLeft(precision.getPrecision());
            do{
                a0=a1;
                b0=b1;
                a1 = a0.add(b0).multiply(oneHalf);
                b1 = bigRoot(a0.multiply(b0),precision);
            }while(a1.subtract(b1).abs().compareTo(TOL)>0);
            context = precision;
            return a1; 
        }

        public BigDecimal getA1(){
            return a1;
        }
        public BigDecimal getB1(){
            return b1;
        }
    }
}
